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ABSTRACT  (MtMimumiOO  words) 

Numerical  solutions  to  the  Navier-Stokes  equations  for  a  3'D,  tiine-dependent,  highly-symmetric 
flow  (Kida,  1985  J.  Phys.  Soc.  Jpn  54)  have  been  completed.  An  effective  resolution  of  up  to 
1024^  collocation  points  (341  modes  after  dealiasing)  is  attained  within  the  memoiy  on  the  256 
MW  CRAY-2  at  Kixtland  AFB  and  the  C}90  at  the  Pittsburgh  Supercomputer  Center.  These 
simulations  constitutes  the  highest  resolution  tuns  made  to  date.  One  of  the  primaiy  purposes  of 
the  work  was  to  create  a  data  base  from  which  a  detailed  energy  transfer  and  triad  analysis  could 
be  made  by  Andrzej  Domaradzki  at  USC  The  data  base  has  l^n  made,  and  runs  for  Reynolds 
numbers  of  500, 1000, 2(XX)  and  5(XX)  have  been  stored  on  tape.  We  shall  give  some  infoimatic 
concerning  the  turbulent  flows  later  in  this  report  The  other  purpose  of  this  work  is  to  try  to 
understand  the  transition  process  through  which  the  flow  becomes  turbulent  Our  early-time 
anaylsis  of  the  data  base  of  runs  was  concerned  with  this  problem,  and  hence  most  of  this  report 
will  deal  with  our  findings.  We  also  attach  a  manuscript  on  this  subject  that  will  be  published 
shortly  in  The  Physics  of  Fluids.  The  transition  process  is  dominated  by  a  thoroughly  unexpected 
and  singular  phenomenon  which  most  certainly  has  ramaflcations  in  the  theory  of  singularities  of 
Euler  and  Navier-Stokes  ^nations,  but  also  may  have  applications  in  turbulence  and  mixing 
away  from  solid  boundaries. 
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Summary 


Numerical  solutions  to  the  Navier-Stokes  equations  for  a  3-D,  time-dependent,  highly-symmetric 
flow  (Kida,  1985  J.  Phys.  Soc.  Jpn  54)  have  been  completed.  An  effective  resolution  of  up  to 
10243  collocation  points  (341  modes  after  dealiasing)  is  attained  within  the  memory  on  the  256 
MW  CRAY-2  at  Kirtland  AFB  and  the  C90  at  the  Pittsburgh  Supercomputer  Center.  These 
simulations  constitutes  the  highest  resolution  runs  made  to  date. 


Goall 

One  of  the  primary  purposes  of  the  work  was  to  create  a  dat?  base  from  which  a  detailed  energy 
transfer  and  triad  analysis  could  be  made  by  Andrzej  Domaradzki  at  USC.  The  data  base  has  been 
made,  and  runs  for  Reynolds  numbers  of  500, 1000, 2000  and  5000  have  been  stored  on  tape.  We 
shall  give  some  information  concerning  the  turbulent  flows  later  in  this  report 


Goal  2 


The  other  purpose  of  this  work  is  to  try  to  understand  the  transition  process  through  which  the 
flow  becomes  turbulent.  Our  early-titne  anaylsis  of  the  data  base  of  runs  was  concerned  with  this 
problem,  and  hence  most  of  this  report  will  deal  with  our  flndings.  We  also  attach  a  manuscript  on 
this  subject  that  is  currently  in  review  in  The  Physics  of  Fluids.  The  transition  process  is 
dominated  by  a  thoroughly  unexpected  and  singular  phenomenon  which  most  certainly  has 
ramafications  in  the  theory  of  singularities  of  Euler  and  Navier-Stokes  equations,  but  also  may 
have  applications  in  turbulence  and  mixing  away  from  solid  boundaries. 
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Background  Information 


□ 


In  order  to  simulate  directly  turbulent  flows  with  high  Reynolds  numbers,  we  make  use  of _ 

symmetries  of  the  Navier-Stokes  equations.  These  symmetries  are  27t  periodicity,  mirror  _ _ 

symmetries  at  planes  xi  =  ntc,  where  i  =  1, 2,  3  and  n  =  0,  ±1,  ±2, ... ,  rotation  (by  ntl)  symmetry  ity  Codes 

Avail  and/or 


Dist 


Special 


around  axes,  (xj.xj)  ®  {nil,  Jt/2)  for  (ij)  =  (1,2),  (2,  3),  (3,  1),  and  a  velocity  permutation 
symmetry  such  that  ui(xi,X2,X3)  =  U2(x2,X3,xi)  =  U3(x3,xi,x2). 

These  symmetries  may  not  be  stable,  and  may  be  observed  in  nature  for  a  short  time  only. 
We  have  found,  however,  that  the  turbulent  flows  that  result  from  these  symmetries  is  qualitatively 
similar  to  flows  which  do  not  enforce  the  symmetries.  Enforcing  the  symmetries,  allows  us  to 
reduce  the  memory  lequiiments  by  almost  200  tioKS  and  the  number  of  operations  per  timestep  by 
almost  1(X).  Simulating  these  prototypical  flows  are  the  only  practical  way  to  study  flows  with  a 
long  range  of  scales  (or  high  Reynolds  numbers).  Otherwise,  we  must  wait  for  over  a  hundred¬ 
fold  increase  in  computing  power  in  order  to  simlate  flows  with  such  a  range  without  the 
symmetries. 

Early  Time  Behavior  and  Transition 

The  transition  to  turbulence  in  this  flow  is  important  for  a  number  of  reasons.  As  in  the  Taylor- 
Green  vortex,  there  is  no  wall  or  forcing  that  allows  a  known  transition  process  to  occur. 
Mechanisms  for  transition  in  these  flows  have  been  conjectiued  to  be  by  vortex  reconnection  or  by 
v(n*tBx-sheet  rollup,  but  have  not  been  proved.  There  are  practical  applications  in  which  a  turbulent 
field  is  desired  in  an  isolated  region  away  from  walls,  say  to  enhance  mixing,  where  a  laminar  one 
exists.  We  have  therefore  tried  to  understand  the  transition  process  in  this  flow. 

The  first  result  is  that  the  transition  is  a  nearly  singular  event  in  space  and  time,  requiring  a  very 
high  resolution.  In  the  above  mentioned  reference  (Kida),  simulations  with  Reynolds  numbers  lA' 
of  up  to  5(XX)  were  presented.  In  our  work,  it  is  shown  that  with  much  higher  resolution,  the 
transition  is  not  well  resolved  and  that  numerical  errors  are  involved  in  the  transition  event  and 
subsequent  evolution.  It  is  interesting  to  note  that  the  transition  event  requires  much  more 
resolution  than  the  turbulence  that  follows.  In  order  to  resolve  properly  and  understand  this 
transition  process,  we  ran  a  succession  of  Reynolds  numbers  from  5(XX)  down  to  500.  A  short 
discussion  of  the  results  are  to  follow;  see  the  manuscript  for  more  details. 

From  simple  low-mode  initial  conditions  and  no  forcing,  the  flow  develops  a  structiue  we  call  a 
"vortex  dodecapole"  around  the  origin.  Twelve  vortex  tubes,  grouped  as  six  dipoles  each 
straddling  an  axis  collapse  symmetrically  towards  the  origin.  This  structure  is  caused  by  the  mirror 
symmetries  and  velocity  permutations  imposed  on  the  flow,  but  may  exist,  for  smne  time,  without 
symmetries.  A  6-vortex  structure  appears  around  the  point  {n/2,.nl2,n/2).  For  some  time  period  in 


the  evolution,  the  dodectqmle  exhibits  what  we  believe  is  a  self-similar  behavior,  collapsing  toward 
the  migin  at  the  same  rate  that  the  tube  radius  is  decreasing  due  m  stretching. 

The  collapse  is  eventually  altered  and  then  halted.  Even  in  the  low  Reynolds  number  case,  the 
highest  resolution  (1024^)  was  required  for  adequate  resolution  (at  least  10  decades  in  range  of 
energy  spectra).  The  energy  spectrum  is  full  but  seems  to  lack  any  particular  power  law  scaling.  In 
the  higher  Reynolds  number  cases,  the  flow  undergoes  transition  in  which  small-scale  vorticity  is 
produced  from  around  the  collapse  points  and  eventually  fills  the  domain  in  an  apparently  random 
fashion.  While  we  believe  that  in  the  high  Reynolds  number  cases,  lack  of  resolution  which  brings 
about  emnrs  and  aliased  modes,  causes  the  dodecapole  to  break  down,  in  careful  observatitm  of  the 
lower  Re  cases,  we  see  what  we  believe  to  be  the  start  of  vortex  reconnection  between  the  tubes  of 
the  dodecapole.  When  partial  reconnection  begins,  the  self  similarity  is  then  lost  and  the  stracture 
becomes  unstable  and  breaks  up  into  small  scale  vorticity. 

While  some  integral  quantities  such  as  energy  and  enstrophy  only  very  weakly  reflect  this  event, 
skewness  and  flatness  show  large  fluctuations.  Local  quantities  like  the  maximum  of  vorticity 
show  even  wilder  behavior.  The  peak  in  enstrophy  occurs  just  before  the  transition.  Skewness 
peaks  sharply  at  the  time  of  the  transition  with  values  of  0.55  fw  Rc=500, 0.70  for  Re=10(X),  0.95 
for  Re«2000,  and  1.1  for  Res50(X).  Skewness  then  drops  to  between  .4  and  .5  in  the  turbulence 
flow  regime.  Flatness  peaks  at  19, 37, 82  and  170,  but  shows  signs  of  nonconvergence.  The  case 
Re=1000  shows  a  difference  in  peak  Harness  with  resolution:  27  for  256^,  34  for  512^,  and  37  for 
1024^.  If  the  peak  values  of  quantities  differ,  the  subsequent  evolution  differs  also  but  to  a  lesser 
extent.  The  maximum  of  vorticity  peaks  spatially  in  the  middle  of  a  vortex  tube  of  the  dodecapole 
and  temporally  at  transition.  Even  in  the  Re=500  case,  the  peak  values  are  100  fra*  256^,  115  for 
512^  and  120  for  1024^.  After  the  breakup  of  the  coherent  vortices,  there  is  a  sharp  drop  in  all 
quantities. 

In  the  higher  Reynolds  number  cases,  it  is  not  clear  whether  viscosity  or  lack  of  numerical 
resolution  is  causing  the  transition  since  similar  but  less  dramatic  behavim-  is  seen  at  lower  Re.  We 
speculate  that  the  dodecapole  may  be  a  finite-time  singularity  of  the  Euler  equations.  The  structure 
of  the  singularity  is  characterized  by  an  analyticity  strip  analysis  and  by  curve  fitting  of  the 
maximum  value  of  vorticity.  The  latter  analysis  shows  that  the  maximum  vorticity  scales  as  (T  - 
Ter)**  shortly  before  the  transition  The  width  of  the  analyticity  strip  (distance  of  the  nearest 
singularity  in  complex  space-time  to  real  space-time)  approaches  zero  faster  than  exponential. 
Exponential  behavior  suggests  no  finite  time  singularity. 


It  is  interesting  to  specul&te  about  experinientally  creating  similar  structures.  When  two  vortex 
rings  are  shot  symmetrically  towards  each  other,  rapid  stretching  occurs  and  after  the  symmetry  is 
broken,  there  is  a  rapid  evolution  of  small  scale  vonicity  and  mrbulence.  The  dodecapole  is 
perhaps  an  extreme  example  of  such  a  phenomenon. 

This  work  has  spawned  other  research  including  my  current  project  on  direct  simulations  of 
turbulence  in  complex  time  in  which  we  attempt  to  establish  the  location  and  form  complex  time 
singularities  in  order  to  find  the  leading  behavior  of  turbulence  in  real  space-time. 

The  Turbulence  Stage 

After  the  peak  in  dissipation  rate  has  been  reached,  the  energy  spectra  begin  to  show  Kolmogorov 
scaling.  We  refer  to  figure  4c  in  the  manuscript  and  note  that  the  k'5/3  scaling  exists  from  a 
wavenumber  of  10  to  150  (1.2  decades).  We  also  refer  to  figure  1 1,  an  isovorticity  surface  plot 
which  shows  an  apparently  random  distribution  of  vorticity  which  gives  rise  to  such  a  scaling. 
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abstract 


The  three-dimensional  time  evolution  of  a  high-symmetry  initial  condition^  is  simu¬ 
lated  using  a  Fourier  pseudospectral  method  for  Re=  Ijv  =  500,1000,2000,  5000  with 
an  effective  resolution  of  1024^  collocation  points  (171^  independent  modes,  maximum 
wavenumber  kmax  —  340).  It  is  found  that  much  before  the  peak  enstrophy  is  reached, 
there  is  a  short  interval  when  the  local  quantities  increase  sharply.  It  is  also  found  that 
during  this  interval,  6  vortex  dipoles  (at  the  origin)  and  3  dipoles  (at  the  7r/2  corner) 
collapse  towards  two  separate  vorticity  null  points  at  the  opposite  corners  of  the  domain 
in  a  .nearly  self-similar  fashion.  The  coherent  vortices  break  up  zdterwards  followed  by 
sharp  decrease  in  local  quantities.  The  singularity  analysis  shows  that,  within  the  limits 
of  the  resolution,  the  maximum  vorticity  scales  approximately  as  (T  —  Tc)”^  shortly 
before  the  break  up.  However,  the  increase  in  peak  vorticity  stops  at  a  certain  time 
possibly  due  to  viscous  dissipation  effects.  The  temporal  evolution  of  the  width  of  the 
analyticity  strip  shows  that  6  approaches  zero  at  a  rate  faster  than  exponential,  but 
reaches  a  minimum  value  and  starts  to  increase.  This  suggests  that  the  solution  remains 
uniformly  analytic  as  is  the  case  in  the  viscous  Burgers  equation. 


1  Introduction 


r- 
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It  is  well-known  that  ^  the  degrees  of  freedom  required  to  simulate  a  three-dimensional 
turbulent  flow  scale  with  and  hence  it  is  extremely  difficult  to  simulate  turbulent 

flows  of  high  Reynolds  number.  The  difficulty  can  be  partially  compensated  by  using 
initial  conditions  having  symmetries.  For  such  flows,  it  is  possible  to  calculate  the 
evolution  in  a  fraction  of  the  whole  flow  field  and  thus  to  increase  the  effective  resolution. 
Navier-Stokes  equations,  in  the  absence  of  solid  botmdaries,  will  preserve  the  symmetries 
of  the  initial  conditioxis  with  time  since  they  are  invariant  under  several  transformation 
groups,  namely,  space  translations,  rotations  and  plane  reflections. 

A  popular  example  of  an  initial  condition  with  high-symmetries  is  the  Taylor- Green 
vortex  By  making  use  of  symmetries,  it  is  possible  to  gain  a  factor  of  4  in  the 
separation  of  scales,  64  in  memory,  32  in  number  of  operations  *.  An  initial  condition 
with  an  even  larger  number  of  symmetries  is  the  high-symmetry  flow  suggested  by  Kida 
The  symmetries  of  this  flow  reduce  the  memory  requirement  by  a  factor  of  192  and  the 
number  of  operations  by  a  factor  of  96  (3  times  more  than  the  Taylor-Green  vortex). 

Both  flows  have  been  analyzed  extensively  by  direct  numerical  simulations.  The 
former  by  Orszag  ®,  Brachet  et  al.  and  the  latter  by  Kida  and  Murakami  *’®.  Table 
1  presents  an  outline  of  the  parameters  of  study  in  these  works.  In  Table  1,  the  last 
three  columns  give  the  resolution  of  the  whole  flow  field  (periodicity  box),  the  first 
one  the  number  of  collocation  points  before  desdiasing  and  the  second  one  the  number 
of  independent  Fourier  modes  after  dealiasing  and  the  last  one  the  maximum  resolved 
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wavenumber  kmax  in  the  flow  field.  To  the  authors’  knowledge,  the  present  work  has  the 
highest  resolution  obtained  in  literature  using  Fourier  spectral  methods.  Also  given  in 
Table  1  is  the  Reynolds  uumber  based  on  the  Taylor  microscale,  Ra,  which  evaluated  at 
the  time  when  the  peak  dissipation  is  attained  (See  Section  2  for  the  definition  of  Ra). 

The  previous  simulations  in  literature  show  that  several  features  of  fully  developed 
turbulence  are  observed  at  late  times  (after  the  peak  in  dissipation  is  reached).  Taylor- 
Green  simulations  show  that  the  energy  spectra  follow  a  power  law  in  wave-number,  the 
exponent  being  between  1.47  to  2.02.  The  Kolmogorov  constant,  on  the  other  hand,  is 
approximately  equal  to  4,  a  value  considerably  larger  than  the  experimentally  observed 
values  of  1.4- 1.5.  Simulations  by  Kida  and  Murakami  give  the  Kolmogorov’s  constant 
to  be  between  1.2  to  2.0,  a  closer  value  to  the  experiments.  Their  energy  spectra  power 
law  exponent  is  calculated  to  be  between  1.5  to  1.7. 

Perhaps  the  more  interesting  stage  in  the  flow  is  the  early  time  evolution;  an  interval 
before  the  peak  in  enstrophy  and  dissipation  is  attained.  Viscous  Taylor-Green  flow 
simulations  show  that  it  is  during  this  interval  the  coherent,  highly-symmetric  vortices 
start  to  break  up.  The  transition  mechanism  is  not  well-understood,  but  there  are  several 
scenarios  proposed  *  such  as  Kelvin- Helmholtz  type  inviscid  instability  mechanisms  or 
viscosity  induced  instabilities  and  possible  vortex  reconnections.  For  the  inviscid  flow 
simulations,  on  the  other  hand,  the  question  as  to  whether  the  smooth  initial  conditions 
will  remain  reguliir  or  not  when  the  three-dimensional  Euler  equations  are  solved  is  yet 
to  be  answered.  Predictions  using  time  series  analysis  and  Pade  approximants  gave 
hints  that  a  singularity  might  develop  spontaneously  in  the  Taylor-Green  flow  and  the 
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structure  of  the  flow  near  the  singular  time  Tc  will  have  important  consequences  for 
the  generation  of  small  scale  turbulent  structures  and  intermittency.  However,  more 
recent  work  by  Brachet  et  al.  ®  using  a  more  extended  series  could  not  confirm  the 
existence  of  a  finite  time  singularity.  An  alternate  approach  "  is  to  examine  the  time 
evolution  of  the  width  of  the  analyticity  strip  S  (the  distance  to  the  real  domain  of 
the  nearest  singularity).  If  6  goes  to  zero  in  finite  time,  the  solution  becomes  singular. 
Simulations  by  Brachet  et  al.  ®  showed  that  6  decreases  exponentially  within  the  limits 
of  resolution  which  suggests  that  the  singularity  will  never  occur  in  finite  time.  However, 
the  resolution  problems  arise  at  an  early  time  leaving  the  results  inconclusive. 

The  present  work  focuses  on  the  high-symmetry  flow  by  Kida  and  Murakami 
particularly  at  the  early  stage  of  the  flow  when  the  coherent  vortices  start  to  break 
up.  Since  the  evolution  of  the  flow  field  in  physical  space  has  not  been  presented  by 
Kida  and  Murakami  *’®,  we  emphasize  the  physical  pictures  (vorticity  contour  sections 
and  isosurfaces)  in  detail.  Section  2  presents  the  initial  conditions,  governing  equations 
and  the  par£uneters  of  study.  Section  3  gives  the  description  of  the  flow  field  with 
time  emphasizing  the  role  of  symmetries  on  the  flow  dynamics.  Section  4  focuses  on 
the  evolution  of  local  and  totzJ  quantities,  with  a  discussion  on  the  possible  singularity 
formation  and  the  role  of  resolution  on  our  results.  Section  5  presents  the  conclusions. 
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2  Equations  Solved,  Initial  Conditions,  Parameters 
of  Study 

The  equations  solved  are  the  three-dimensional  Navier-Stokes  equations  in  rotational 
form  given  by 

^  —  w  X  w  =  —  vn -I- (1) 

with  the  incompressibility  condition  given  by 

V  •  u  =  0  ,  (2) 

where  u  and  u  are  the  velocity  and  vorticity  vectors  respectively  and  11  =  ^  -j-  ( 1  /2)ti^ 
and  the  density,  /)  =  1  , 

The  initial  conditions  are  given  by 

u{x,  y,2,t  =  0)  =  sin  x(cos  3y  cos  z  —  cosy  cos  3z)  (3) 

v(x,  y^z,t  =  0)  =  sin  y{cos  Sz  cos  x  —  cos  z  cos  3x)  (4) 

w{x,  y,z,t  =  0)  =  sin  2(cos  3x  cos  y  —  cos  x  cos  3y)  (5) 

where  u,  v,  w  denote  the  velocity  vector  components  in  x,  y,  z  directions.  This  flow  has 
been  suggested  by  Kida  and  Murakami  and  a  detailed  summary  of  the  flow  symmetries 
Me  discussed  by  these  authors  In  this  work,  we  will  refer  to  the  initial  condition  as 
the  ‘high-symmetry  initial  condition’. 

The  symmetries,  which  we  shall  only  list  here,  occur  in  all  three  directions,  and  are 
a  2?r  periodicity,  a  bilateral  symmetry  through  pleines  htt,  n  =  0,  ±1,  ±2, . . .,  a  7r/2 
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rotational  symmetry  around  the  axis  (7r/2,  t/2),  and  a  permutation  symmetry  of  the 
velocity  components,  u{x,y,z)  =  v{y,z,x)  =  w{z,x,y). 

If  the  above  flow  is  compared  to  the  Taylor-Green  initial  conditions  in  Fourier  space, 
certain  differences  are  found.  The  number  of  non-zero  Fourier  coefficients  of  the  velocities 
for  the  Taylor-Green  flow  is  2,  compared  to  6  for  the  high-symmetry  initial  condition. 
The  magnitude  of  the  initial  wave-number  vector  k  corresponding  to  non-zero  Fourier 
modes  of  the  velocities  is  -|- 1  -t- 1  =  for  the  Taylor-Green  flow  compared  to  ■v/IT 
for  the  high-symmetry  initial  condition. 

Figures  1  (top  left)  and  2  (top  left)  give  an  isosurface  of  the  magnitude  of  the  vorticity 
u  =  O.ScJo  of  the  initial  condition  from  two  views. 

We  solved  the  governing  equations  using  a  Fourier  pseudo-spectral  method  with 
the  “2/3  rule”  dealiasing  We  used  the  second-order  accurate  method  of  time-split 
fractional  steps  for  time  integration.  The  details  of  the  numerical  method  can  be  found 
in  our  previous  work  Appendix  B. 

The  even-mode,  odd-mode  cosine  and  sine  tr£insforms  are  evaluated  using  an  FFT  for 

I 

j  read  sequences  and  suitable  pre-  and  post-processing  as  given  by  Cooley  et  al.  and 

I 

Brachet  et  al.  ®.  Only  four  arrays  and  four  three-dimensional  transforms  are  required 
per  time  step.  An  in-place  indicial  permutation  to  obtain  the  three  velocity  components, 
j  when  needed,  is  used. 

I  The  parameters  of  the  study  are  presented  in  Table  2.  It  is  worth  emphasizing  what 

i  is  meant  by  resolution  in  this  table.  As  an  example,  for  Run  B2,  the  total  resolution  is 

i  given  by  1024^  which  is  the  total  number  of  collocation  points  in  the  periodicity  box. 
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This  is  also  equivalent  to  256  independent  Fourier  inodes.  We  use  dealiasing  in  our  code 
which  reduces  the  independent  Fourier  modes  to  2/3  x  256  =  171  (and  a  maximum 
wavenumber  kmax  =  341).  Note  that  the  symmetries  allow  us  to  solve  the  problem  in 
the  fundamental  box  which  has  a  resolution  of  (1024/4)^  =  256^. 

As  seen  from  Table  2,  we  have  four  runs  with  Reynolds  numbers,  {Re  =  l/v)  500  for 
Rim  C,  1000  for  Run  D,  2000  for  Run  A  and  5000  for  Run  B.  For  Run  A,  remeshing  is 
done  at  t=0.785  by  appending  zero  Fourier  modes,  thus  increasing  the  total  resolution 
from  256^  to  512^,  and  Run  A2  is  started  at  this  time.  Further  remeshing  is  done  on 
Run  A2  at  t=1.625  increasing  the  resolution  from  512^  to  1024^  and  Run  A3  is  started 
at  this  instant.  We  also  continued  both  Runs  A1  and  A2  until  t=2.25  with  the  lower 
resolution. 

For  Runs  C  and  D,  remeshing  is  done  at  t=0.8  increasing  the  total  resolution  from 
256®  to  512®,  and  Run  C2  (and  D2)  is  started  at  this  time.  Further  remeshing  is  done 
on  Run  C2  (and  D2)  at  t=1.55  increasing  the  resolution  from  512®  to  1024®  and  Run 
C2  (and  D2)  is  started  at  this  time. 

For  Run  B,  we  started  with  a  higher  resolution,  512®  and  increased  the  resolution  to 
1024®  at  t=1.625  and  started  the  Run  B2. 

Later,  it  will  be  clear  to  the  reader  after  reading  Section  III  and  IV  why  we  stop  the 
nms  at  the  particular  times.  We  postpone  this  discussion  to  these  sections. 

For  each  resolution,  the  time  step  is  determined  from  the  CFL  condition  and  checked 
continuously  throughout  each  run.  It  was  possible  to  use  an  even  larger  At  for  Runs  Al, 
Cl,  D1  but  we  kept  it  the  same  as  in  Runs  A2,  C2,  D2  since  the  computational  cost  of 
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Runs  Al,  Cl,  Dl  was  very  low. 


The  computational  time  and  memory  requirements  of  the  runs  are  presented  in  Table 
III.  We  used  the  (single  processor)  Cray  2  at  the  Kirtland  Supercomputer  Center  and  the 
Cray  Y-MP  at  the  Pittsburgh  Supercomputing  Center  for  one  low  resolution  run.  Cray 
Y-MP  is  more  than  twice  faster  than  the  Cray  2  but  does  not  have  sufficient  memory  for 
the  high  resolution  runs.  Also,  the  turnaround  time  at  the  Kirtland  Center  was  much 
faster  than  in  Pittsburgh.  For  some  of  the  post-processing  calculations,  and  short  runs, 
we  used  the  newly  installed  Cray  90  (C'90)  at  the  Pittsburgh  Supercomputing  Center 

We  calculated  certain  total  and  local  quantities  in  each  run  during  the  evolution. 
These  were  stored  at  every  t=0.0125.  We  stored  the  velocity,  vorticity  and  dissipation 
vector  fields  at  every  t=0.125  and  sometimes  every  t=0.0625. 

We  calctilated  the  total  energy  E  ;  and  total  enstrophy  ft  ;  given  by 


1 

II 

f  +  u*  +  w^)dV 

(6) 

ft  =  —  i 
V  J 

(7) 

where  Wx,  Wy,  u;*  are  the  vorticity  vector  components  in  x,  y,  z  directions  respectively. 
For  the  discrete  problem,  the  integrals  are  sununations  over  all  grid  points  in  the  solution 
domain  and  V  is  the  total  number  of  grid  points. 

Also  calculated  are  the  total  dissipation  e,  the  velocity  derivative  skewness  S  2md 
flatness  F  given  by 

€(<)  =  (8) 

_  <  {dujdxf  > 

<  {dnidxf  >3/2 
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_  <  {dujdx)*  > 

<  {dujdxY  >2 

where  <  >  denotes  averaging  over  all  grid  points  at  a  given  time,  and  S,j  is  the  strain 

rate  tensor.  We  emphasize  that  our  skewness  and  flatness  calculations  are  performed 
directly  using  the  above  relations  without  making  any  additional  assumption  as  isotropy. 


We  also  calculated  the  Reynolds  number  based  on  Taylor  microsc£ile  Ra,  Kolmogorov 
microscale  17  and  Kolmogorov  wavenumber  kr,  defined  as 


0 

II 

(11) 

(12) 

K  - 1/1 

(13) 

Several  other  local  quantities  are  calculated.  These  are  the  local  maximum  velocity, 
vorticity,  dissipation  Umo*,  Cmox* 

And  finally,  we  calculated  the  band-averaged  three  dimensional  energy  spectrum 
E(k,t)  given  by 

=  (») 

^  ktA 

where  A  is  the  domain  for  p  —  A  <  fc  <  p  -1-  A,  p  =  0, 1, . Veilues  for  A  are  1  /2  and  1. 

We  denote  the  magnitude  of  the  wavenumber  vector  by  k,  that  is;  k  =  |J:|. 


3  Flow  Evolution 

We  start  this  section  giving  the  statistics  of  Run  Bl.  This  run  has  the  same  initial 
conditions  as  one  of  the  rims  Kida  and  Murakami  ^  simulated  so  we  will  compare  our 
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results  with  that  run  in  detail. 


3.1  Late  Stages 

The  evolution  of  energy  E,  enstrophy  fl,  dissipation  c  and  the  Taylor  Reynolds  number 
R\  are  given  in  Figure  3.  We  compare  these  results  with  Figures  1,  2,  3  and  4  of  the 
work  by  Kida  and  Murakami  We  find  that  there  is  perfect  agreement  between  the 
two  simulations. 

As  seen  from  Figure  3,  the  energy  stays  constemt  until  t=2  showing  that  the  flow  is 
close  to  inviscid  at  this  stage.  The  enstrophy  and  dissipation  increase  smoothly  (without 
any  §udden  jump),  the  log-linear  plot  shows  an  almost  straight  line  with  time  implying 
that  the  initial  enstrophy  increase  is  at  least  exponential.  In  fact,  there  a--  indications 
that  the  almost  straight  line  in  Figure  3  for  enstrophy  is  slightly  parabolic  which  would 
give  a  faster  than  exponential  growth  rate.  A  final  observation  in  Figure  3  is  that  there 
is  a  small  dimple  observed  in  E,  and  Rx  between  t=2  and  3.  This  clearly  shows  a 
change  in  the  rate  of  increase  (decrease)  but  at  this  stage  it  is  not  clear  why  it  is  so. 

The  flow  field  after  the  peak  e  is  attained  at  t=3.5  starts  showing  statistics  in  agree¬ 
ment  with  Kolmogorov’s  scalings.  Figure  4a  and  4b  shows  the  spectra  eind  the  spectra 
multiplied  by  A:®/®  at  t=6.  We  observe  Kolmogorov  scaling  between  k=10  to  70.  There 
is  a  very  slight  bump  at  k=170  which  is  noticeable  but  not  severe  enough  to  suggest  the 
resolution  problems.  Figure  4c  shows  the  spectra  multiplied  by  A:®/®  for  Run  B2  (our 
highest  resolution  run)  for  different  late  times  (6  different  times  between  t=5.75  and 
6.00)  which  exhibits  Kolmogorov  scaling  between  the  wavenumbers  k=10  to  150. 
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Questions  arise  as  to  how  the  coherent  initial  vortices  evolve  into  isotropic  structures 
giving  Kolmogorov  scaling,  what  the  transition  mechanisms  are,  when  the  transition 
takes  plaice  and  whether  there  is  enough  resolution  in  this  high  Reynolds  number  run  or 
not.  One  resolution  adequacy  check  is  to  calculate  the  time  evolution  of  the  Kolmogorov 
wavenumber  kr,  amd  compare  it  to  the  highest  resolved  mode  in  Run  B1  (which  is  170). 
Figure  5  gives  this  evolution.  We  see  that  k„  increases  to  323,  which  is  above  the 
highest  resolved  mode  in  our  run.  We  conclude  that  resolution  has  to  be  increased 
somewhere  between  t=l  and  2  and  both  our  Run  B1  and  Kida-Muraikami  runs  axe 
probably  underresolved.  This  is  the  maun  reason  of  our  doubling  the  resolution  at 
t=1.625  and  starting  Run  B2  which  has  a  maximum  wavenumber  of  341. 

To  investigate  the  transition  mechanisms  to  turbulence,  we  return  to  the  flow  pictures 
in  physical  space  in  the  next  section. 

3.2  Early  Stages 

3.2.1  Symmetries  and  Dynaunics 

In  order  to  understand  the  eaxly  time  evolution  of  the  flow,  we  focus  on  two  sepairate 
regions  in  the  fundaunental  box  where  symmetries  result  in  the  formation  of  interesting 
vortex  structures.  These  regions  au:e;  i)  x,  y,  z  =  x/2  faces,  ii)  x,  y,  z  =  0  faces.  We 
will  focus  on  the  XY  plaines  at  z=x/2  and  z=0  to  describe  these  symmetries.  Similau’ 
conclusions  cam  be  drawn  for  the  other  two  plames. 
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For  the  XY  plane  z=jr/2  the  following  symmetries  exist: 


u(x,y,7r/2)  =  v(y,x,ir/2) 

(15) 

u?(i,y,jr/2)  =  -t£>(y,  X,  5r/2) 

(16) 

A  similar  symmetry  is  also  present  in  the  vorticity  field. 

This  is  not  written  explicitly 

by  Kida  and  Murakami  but  can  be  easily  seen  by  taking  the  curl  of  the  initial  velocity 

field.  Then,  it  can  be  shown  that: 

y,  Jr/2)  =  x,  z/2) 

(17) 

uJz(x,y,ir/2)  -  -a;,(y,  x,  t/2) 

(18) 

It  is  seen  from  the  above  equation  that  the  normal  vorticity  (a;,  in  plane  XY)  is  an 
odd  ftmction  with  respect  to  x=:y.  This  is  nothing  but  a  necessary  condition  to  obtain 
a  dipole  along  the  line  x=y.  The  dipole  is  not  necessarily  aligned  with  the  z  axis,  in 
addition  to  its  out-of-plane  vorticity  u;,  it  also  has  non-zero  vorticity  components  in  x 
and  y  directions. 

Before  giving  the  picture  of  the  dipole,  we  emphasize  its  one  final  feature:  At  t=0, 
Ux{x,y,irf2)  =  — u>*(y,z,7r/2)  =  0,  thus  there  is  no  out-of-plane  vorticity  in  the  initial 
field.  However,  there  is  no  restriction  that  the  vorticity  will  stay  zero  at  later  times  at 
z  =  7r/2  plane.  In  fact  we  see  in  Figure  6  that  this  is  the  case. 

Similar  sjrmmetries  exist  for  the  other  planes.  There  is  another  dipole  in  the  y  =  7r/2 
plane  along  x=z,  and  still  another  in  the  plane  x  =  z/2  along  y=z.  These  three  dipoles 
are  oriented  in  such  a  way  that  they  all  move  towards  the  point  x,  y,  z  =  ir/2, 7r/2, 7r/2. 

The  other  regions  of  interest  are  the  x,  y,  z  =0  planes.  We  will  consider  z=0  plane 
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to  describe  the  symmetries.  For  this  plane,  we  have 


w*(x,y,0)  =  -a?*(-x,y,0) 

Wj:  =  0 
Wj,  =  0 


(19) 

(20) 
(21) 


We  see  that  the  out-of-plane  vorticity  is  an  odd  function  with  respect  to  the  y  axis 
and  thus  its  image  with  respect  to  the  y  axis  will  form  a  dipole.  Similarly,  we  have 
another  symmetry  with  respect  to  the  x  axis.  This  is  given  by 


(22) 


Then,  due  to  the  above  odd  symmetry  with  respect  to  the  x  axis,  the  image  of  the  vortex 
with  respect  to  the  x  axis  will  form  another  dipole. 

The  symmetries  have  an  additional  property  at  the  zero  faces;  The  normal  vorticity 
is  zero  along  the  zero  edges,  and  is  zero  at  the  origin.  Then,  the  origin  is  a  vorticity  null 
point. 

With  all  these  symmetries  in  mind,  we  describe  the  evolution  of  the  high-symmetry 
flow  in  the  next  section. 


3.2.2  Evolution 

Figures  1,  2,  7-12  display  magnitude  of  vorticity  isosurfaces  at  different  times  for  Run 
B  £rom  two  different  angles.  The  first  view  angle  is  used  in  Figures  1,  7,  9,  11,  12  and 
shows  the  x,y,z  =  ir/2  planes  clearly.  The  second  view  angle  is  used  in  Figures  2,  8,  10 
and  shows  the  origin. 
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The  initiad  stages  of  the  evolution  is  seen  in  Figures  1,  2.  Figure  1  shows  that  three 
dipoles  are  immediately  formed  with  oppositely-signed  vortices  appro<iching  towards  one 
another  due  to  the  strain  field  set  up  in  the  flow  field.  At  t=0.75  (Figure  1,  bottom  right), 
the  standard  head-tail  cross  sections  of  the  three  dipoles  in  the  x,y,2  =  ‘k/2  pl2mes  are 
clearly  seen.  In  the  meantime,  a  different  motion  takes  place  near  the  origin.  To  see 
the  dynamics,  we  consider  the  smaller  vortex  given  in  Figure  2  (bottom).  This  vortex 
will  have  an  image  vortex  located  below  with  opposite  sign  (quadrant  4),  another  image 
vortex  located  to  its  left  of  opposite  sign  (quadrant  2)  and  finally  the  image  vortex  in 
quadrant  2  will  have  an  image  in  quadrant  3.  Therefore,  totally,  we  have  4  vortices  or 
2  dipoles  in  the  z=0  plane.  A  similar  argument  can  be  made  for  x=0  and  y=0  planes 
giving  8  more  vortices  or  4  more  dipoles.  So,  there  are  12  vortices  or  6  dipoles  on  the 
zero  faces  approaching  the  origin.  We  shall  call  this  configuration  a  dodecapole  since  we 
can  treat  it  as  an  entity  (six  dipoles)  straddling  each  of  the  axes.  We  emphasize  that 
in  all  vorticity  isosurfaces  given  in  this  work,  only  the  fundamental  box  is  displayed. 
Therefore,  only  three  vortices  (1/4  th  of  the  12)  will  be  seen  in  the  pictures. 

The  three  vortices  of  the  dodecapole  approach  to  the  edges  and  start  to  get  flattened. 
The  flattening  can  be  better  seen  in  Figure  13  which  shows  the  magnitude  of  vorticity 
contours  at  t=0.75, 1.0, 1.25  and  1.5  for  Run  A.  The  topology  looks  similar  to  the  typical 
flattening  and  subsequent  head-tail  formation  of  a  dipole.  (For  example:  See  Figures 
6  and  7  by  Kerr  '^).  We  remind  the  reader  that  our  configuration  has  three  flattened 
vortices  in  each  plane  (or  12  when  the  images  are  considered)  converging  to  a  null  point 
(the  origin). 
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Figures  7  and  8  show  the  vorticity  isosurfaces  at  t=1.25,  1.5,  1.75,  1.81.  The  vorticity 
has  amplified  considerably  at  this  stage  and  the  threshold  is  raised  for  clarity.  Three 
dipoles  directed  towards  the  point  jr/2,  ?r/2,  7r/2  is  clearly  seen  at  t=i.25  (Figure  7, 
top  left).  At  t=1.5,  more  vortex  tubes  appear  on  each  ir/2  face  (Figure  7,  top  right). 
When  part  of  each  main  dipole  exits  on  a  face,  it  enters  on  two  other  faces  due  to  the 
symmetries  which  are  given  by 

w*(jr/2,  y,  z)  =  -(Ajy(x,  x-/2,  z)  =  -a;,(r,  y,  ff/2)  (23) 

At  t=1.75  and  1.81  (Figure  7,  bottom),  the  dipole  structure  is  becoming  distorted  and 
the  flow  is  becoming  increasingly  complicated.  To  see  the  complexity,  we  display  in 
Figure  14,  a  contour  section  at  z  =  x-/2  for  the  times  t=1.25,  1.5,  1.625,  1.75.  (Recall 
the  symmetry  presented  before  with  respect  to  x=y  line.)  In  the  meantime,  the  vortices 
of  the  dodecapole  continue  to  get  flattened  as  seen  from  Figure  8.  The  threshold  chosen 
does  not  show  the  connectivity  between  the  dodecapole  and  the  dipoles  but  the  main 
reason  of  the  choice  is  to  focus  on  the  stn^cture  near  the  origin. 

Figures  9  and  10  display  the  evolution  at  t=1.875,  1.937,  2.0  and  2.0625.  The 
'  vorticity  threshold  is  raised  even  more  compared  to  the  previous  figures.  It  seems  that 
the  vorticity  is  becoming  smaller  in  scale.  Both  the  dipoles  approaching  the  r/2,  fl'/2, 
x/2  comer  and  the  three  vortices  of  the  dodecapole  approaching  the  origin  start  breaking 
up.  A  vorticity  cross  section  given  in  Figure  15  shows  that  the  structure  still  preserves 
its  shapes  until  t=2.062  but  afterwards  (Figure  16),  tremendous  sunount  of  small-scale 
vorticity  is  observed. 

From  all  these  pictures,  it  is  clear  that  a  major  break  up  t^lkes  place  for  Run  B 
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roughly  between  t=1.875  and  t=2.5.  We  have  not  displayed  the  isosurfaces  for  runs  A, 
C  and  D  but  these  runs  also  have  qualitatively  simileur  pictures  in  the  same  interval. 
There  are  slight  quantitative  differences  which  are  discussed  in  Section  4.1. 

Now  if  we  return  to  Figure  3,  we  observe  that  none  of  the  quantities  displayed  show 
any  major  change  during  this  interval.  Perhaps  enstrophy  and  dissipation  rate  changes 
slightly  but  this  is  hardly  noticeable.  It  is  well-known  that  most  of  the  integrated 
quantities  in  vortex  interactions  are  not  sensitive  enough  to  reflect  sudden  changes  in 
the  flow  field  due  to  the  violent  intermittent  events.  However,  there  are  some  exceptions 
such  as  velocity  derivative  skewness  and  flatness  which  reflect  the  sudden  changes  in  the 
vortex  interaction  dynamics  properly.  For  example,  in  vortex  reconnection  problems, 
when  the  interaction  is  at  its  most  intense  stage,  sharp  increases  in  skewness  values  are 
observed  A  similar  trend  is  seen  in  the  work  by  Bracket  et  al.  ®  (their  Figiu*e  10)  and 
Domaradzki  et  al.  (their  Figure  Ic). 

In  order  to  quantify  the  transition  dynamics,  we  look  at  the  evolution  of  certain 
qusmtities  which  could  capture  the  important  events.  That,  we  present  in  the  next 
section. 


4  Transition  Stage:  Quantification 

4.1  Evolution  of  Local  and  Total  Quantities 

We  present  in  Figures  17-18  the  evolution  of  total  and  local  quantities  as  energy  E,  en¬ 
strophy  fl,  Kolmogorov  wavenumber  A:,,  Taylor  Reynolds  number  Rx,  velocity  derivative 
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skewness  S  and  flatness  F  and  maximum  vorticity  Umax-  Due  to  the  space  limitation,  we 
will  only  present  the  results  of  the  highest  Reynolds  number  run  (Run  B).  The  results 
obt2uned  in  Runs  A,  C  and  D,  will  be  mentioned  whenever  necessary.  Figures  17  and 
18  give  results  from  Run  B  with  a  comparison  of  the  results  of  two  different  mesh  sizes, 
Run  B1  (dashed)  and  Run  B2  (solid).  Figure  19  gives  a  comparison  of  two  of  the  runs 
having  higher  resolution,  Run  A3  (solid)  and  Run  B2  (dashed). 

We  see  in  Figures  17  and  19  that  evolution  of  energy,  enstrophy,  Kolmogorov  wavenum¬ 
ber  (which  is  also  a  function  of  enstrophy)  and  Taylor  Reynolds  number  (which  is  a 
function  of  enstrophy  and  energy)  axe  not  significantly  different  for  different  resolutions. 
Quantitatively,  if  Run  B2  and  B3  are  compared,  the  energy  differs  by  less  than  1%, 
enstrophy  less  than  5%. 

Next,  we  discuss  the  evolution  of  certain  local  quantities.  We  see  from  Figure  18  that 
approximately  between  t=:1.5  and  2.5,  there  are  drastic  changes  in  the  local  maximum 
vorticity.  In  this  short  interval,  local  maximum  vorticity  uJmax  increases  13  times  for 
Run  A3,  16  times  for  Run  B2  (Figure  18,  bottom),  6.5  times  for  Run  C3  and  10  times 
for  Run  D3.  Skewness  increases  as  much  as  1.1  for  Run  B2  (Figures  18  and  20),  0.9  for 
Rtm  A3,  0.55  for  Run  C3  and  0.70  for  Run  D3  .  Flatness  increases  to  170  for  Run  B2 
(Figures  18  and  20),  80  for  Run  A3,  19  for  Run  C3  and  37  for  Run  D3.  The  smallest 
change  is  found  in  the  local  maximum  velocity  Umax  which  increases  less  than  twice  in 
Run  B2  (not  shown). 

All  these  sudden  increases  are  followed  by  sudden  drops  which  are  as  sharp  as  the 
increase  rates.  Certain  quantities  start  to  show  fluctuations  after  this  stage  (such  as 
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Mmox  wd  S  in  Figure  18). 


Now  it  becomes  clear  that  the  display  of  these  local  quantities  illuminates  a  violent, 
intermittent  event  taking  place  in  the  flow  fleld.  In  the  absence  of  this  diagnostics  (as  is 
the  case  in  Kida  and  Murakami  ^),  it  is  not  possible  to  see  the  intensity  of  this  event  by 
only  looking  at  energy  and  enstrophy  pictures  only.  The  event  is  so  sudden  and  intense 
that  all  the  local  fimctions  appear  to  look  like  the  delta  function  at  this  instant. 

Questions  arise  as  to  how  effectively  this  large  local  qu£intity  can  be  resolved  and  how 
much  the  results  vary  with  resolution.  It  is  evident  that  for  the  extreme  hypothetical  case 
when  one  of  the  local  quantities  is  a  delta  function,  it  can  be  only  resolved  with  infinite 
resolution.  An  ideal  shock  wave  (with  zero  thickness),  for  example,  in  the  absence  of 
any  smoothening,  needs  infinite  number  of  grid  points  to  be  resolved. 

In  appendix  A,  we  show  the  dependence  of  the  Fourier  transform  of  a  delta  function 
on  the  munerical  resolution  (the  cut-off).  It  is  seen  from  this  derivation  that  for  such  a 
function,  doubling  the  resolution  will  double  the  amplitude  of  the  function. 

For  Runs  A  and  B,  we  find  that  the  local  peaks  axe  very  large  and  as  the  resolution  is 
doubled,  the  value  of  the  peak  fotmd  in  the  higher  resolution  run  is  approximately  twice 
as  much  as  in  the  lower  resolution  run.  Unfortunately,  it  is  not  currently  possible  to  go 
above  a  resolution  of  1024®  due  to  severe  limitations  in  computer  time  and  memory.  On 
the  other  hand,  for  Runs  C  and  D,  there  is  a  convergence  in  local  quantities  within  10 
to  20  %  respectively  (between  C2  and  C3,  D2  and  D3).  A  detailed  resolution  analysis 
is  presented  in  Appendix  B. 

Such  a  sudden  increase  in  the  local  quantities  in  the  flow  field  naturally  brings  the 
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singularity  arguments  into  the  picture.  We  investigate  the  issue  in  detail  in  the  next 
section. 


4.2  Singularity  Analysis 

The  regularity  and  uniqueness  theorems  for  the  three  dimensional  Euler  and  Navier- 
Stokes  equations  are  far  from  being  complete  in  literature.  An  important  result  is  by 
Constantin  who  shows  that  the  breakdown  of  smooth  solutions  to  3-D  incompress¬ 
ible  Navier-Stokes  equations  cannot  occur  without  the  corresponding  solution  of  the 
incompressible  3-D  Euler  equations.  Therefore,  a  possible  singularity  formation  in  our 
problem  will  have  important  implications  on  the  formation  of  a  singularity  in  3D  Euler 
equations. 

We  start  with  the  evolution  of  during  the  time  interval  t=1.64-2.25  given  in 
Figure  20.  The  top  figure  shows  the  trend  in  two  runs  (B2  and  A3)  shortly  before  the 
peak  vorticity  is  reached.  We  observe  three  intervals  with  different  slopes  in  this  figure. 
The  first  slope  increase  is  around  t=1.78  for  Run  B2  (t=1.88  for  Run  A3),  £ind  another 
increase  around  t=1.93  for  Run  B2  (t=1.96  for  Run  A3).  If  the  evolution  continues  with 
this  p2ice,  it  is  clear  that  the  function  will  hit  the  time  axis  amd  the  local  vorticity  will 
diverge.  We  performed  a  least  square  fit  in  Figure  20  for  the  regions  with  the  steepest 
slope  and  fotmd  that  the  time  axis  will  be  intersected  at  Tc  =  2.0389  for  Run  B2  and 
Tc  =  2.1077  for  Run  A3.  However,  as  seen  from  Figure  20  (bottom),  the  slope 
gets  less  steep,  reaches  a  local  minimum  and  starts  to  increase.  We  repeated  the  same 
analysis  for  our  lower  Reynolds  number  runs  C  and  D  and  obtained  similar  results  as 
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seen  in  Figure  21.  Extrapolation  of  the  curves  to  the  x  axis  gives  Tc  =  2.21  for  Run 
C3,  and  Tc  =  2.45  for  Run  D3.  Note  that  the  slope  in  Figure  21  becomes  steeper  as 
approached  to  Tc  similar  eis  in  the  runs  A  and  B. 

It  is  important  to  investigate  how  Umax  and  (Tc  —  T)  scale  particularly  during  the 
interval  with  the  steepest  slope  in  Figure  20.  This  is  presented  in  Figure  22  (recall  that 
the  Tc’s  used  are  not  the  same  for  two  runs).  Least  squares  fitting  for  Run  A3  between 
t=1.9640-2.0515  (using  8  data  points)  gives  u  ~  (Tc  —  T)“*  ®®  £ind  for  Run  B2  between 
t=  1.9265-1.9890  (using  6  data  points)  gives  u  ~  (Tc  —  T)"^  “.  The  power  laws  found 
depend  slightly  on  the  number  of  data  points  used  in  the  fitting.  This  is  a  consequence 
of  the  fax:t  that  the  points  are  not  located  exactly  on  a  straight  line  but  slightly  curved. 
Taking  various  intervals  and  changing  the  number  of  points  gave  exponents  between 
0.9983  and  1.181  for  all  runs  but  never  close  to  2.  We  emphasize  the  fact  that  the 
highest  resolution  results  are  used  to  calculate  the  scaling  and  it  is  not  possible  to  go 
above  a  resolution  of  1024^  to  further  check  the  resolution  dependance  of  the  exponents. 

The  scaling  exponent  found  is  similar  to  the  scaling  calculated  by  Kerr  in  his 
single  dipole  simulation.  Recently,  Bhattacharjee  and  Wang  have  suggested  an  ana¬ 
lytical  model  of  three-dimensional  Euler  flows  which  gives  a  singularity  on  a  line  joining 
two  vorticity  nulls.  It  is  found  that  vorticity  diverges  as  (Tc  —  T)”^.  A  similar  scaling  is 
calculated  by  Grauer  and  Sideris  in  axisymmetric  flows  with  swirl.  However,  Pumir 
and  Siggia  find  a  different  scaling  for  a  similar  problem  with  higher  resolution  which 
is  calculated  to  be  (Tc  —  T)~^. 

Next,  we  examine  the  time  evolution  of  the  width  of  the  analyticity  strip  6.  The  width 
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6  is  the  distance  of  a  complex-domain  singularity  to  the  real  axis.  It  ctin  be  determined 
from  the  large  wavenumber  behavior  of  the  energy  spectra.  For  large  wavenumbers,  it 
can  be  shown  that  the  leading  term  of  the  spectra  is  given  by: 

E{k,t)  ~  exp{—2S{t)k)  (24) 

if  6  remains  finite  in  time.  If,  however,  the  width  6  decreases  to  zero  in  finite  time,  a 
complex-domain  singularity  reaches  the  real  axis  and  the  spectra  behaves  like  a  power 
law  in  wavenumber  given  by 

(25) 

We  calculated  the  temporal  evolution  of  S  in  all  our  simulations  by  performing  least 
squares  fitting  assuming  two  different  shape  functions  for  the  spectra.  The  first  (we  will 
refer  to  it  as  Type  1  fitting)  is  in  the  form  as  exp(— 2£(<)fc)  and  the  higher  end  of 
wavenumber  interval  100  <  A:  <  340  is  used  in  the  fittings  for  all  runs.  The  second  (we 
wiU  refer  to  it  as  Type  2)  is  in  the  form  Ak^  exp{—26{i)k)  and  the  fitting  is  done  for 
wavenumber  intervals  3.5  <  k  <  340  and  9.5  <  k  <  340.  The  results  obtained  from  these 
two  intervals  are  almost  the  same  for  delta  values.  This  can  be  seen  in  Figure  23  which 
gives  the  plots  of  the  energy  spectra  of  Run  C3  at  t=2.3  and  two  fit  curves  using  the 
second  type  of  fitting  using  the  intervals  3.5  <  k  <  340  (solid  curve)  and  9.5  <  k  <  340 
(dashed  curve).  As  seen  from  this  figure,  the  fit  curve  follows  the  points  obtained  from 
the  simulation  perfectly. 

The  energy  spectra  for  two  different  times  (t=1.75  and  2.25)  for  Run  B2  are  presented 
in  Figure  24a.  Figure  24b  displays  k^E{k)  and  Figure  24c  displays  k'^E{k)  .  Note  that 
Figure  24c  implies  that  somewhere  between  t=1.75  and  2.25,  the  spectra  will  be  nearly 
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horizontal  line  which  indicates  a  seeding  approximately  as  k~*.  Recall  that  at  this 
instant,  the  whole  flow  field  (fundamental  box)  has  3  dipoles  near  the  origin  and  3  near 
the  7r/2  corner.  The  dipoles  can  be  considered  as  flattened  vortex  sheets  with  sudden 
vorticity  jumps  across  them  (vorticity  changes  sign  across  the  flattened  dipole  sheet). 
The  configuration  is  similar  to  SaiFman’s  vorticity  field  with  random  two-dimension2d 
vorticity  discontinuities  which  also  gives  a  k~*  scaling.  As  the  structure  breaks  down, 
the  spectra  evolves  towards  k~^  but  continues  to  change.  Therefore,  k~^  scaling  might 
be  oidy  a  tr£Uisient  intermediate  stage  as  the  cascade  evolves  into  the  Kolmogorov  regime 
at  later  times  (A.Pouquet,  private  communication). 

Returning  to  the  6  argument,  we  present  the  temporal  behavior  of  the  analyticity 
width  in  Figure  25  (Type  2  fitting),  and  display  the  results  in  log* log  scale.  The  6  values 
found  using  the  Type  1  fitting  are  found  to  be  qualitatively  the  same  and  quantitatively 
slightly  larger  that  those  calculated  by  T3rpe  2  fitting.  The  horizontal  dashed  line  in 
Figure  25  shows  twice  our  effective  resolution. 

One  interesting  observation  from  Figure  25  is  the  existence  of  the  regions  of  different 
slopes  for  all  four  runs.  For  Run  A3,  after  the  first  sharp  drop  in  6,  there  is  a  region 
of  almost  constant  slope  between  t= 1.789- 1.9765  (recall  that  our  predicted  Tc  =  2.1077 
for  this  run)  and  another  region  with  a  steeper  slope  between  t= 1.98-2.1  after  which  the 
curve  hits  a  minimum  value  at  t=2.1265.  The  width  starts  to  increase  after  this  time. 
A  similar  trend  is  seen  for  Run  B2.  After  the  initial  wiggly  region,  the  curve  shows  a 
constant  slope  between  t=l. 7765-1. 8265  (recall  that  the  predicted  Tc  =  2.0389  for  Run 
B2).  After  t=1.83,  it  seems  like  the  decrease  in  6  slows  down.  This  slow  decrease  in 
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6  continues  until  t=1.964.  After  this  time,  there  is  a  decrease  in  the  slope  which  stays 
uniform  between  t= 1.964-2.064.  At  t=2.064,  the  curve  reaches  to  its  minimum  value. 


% 
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Similuly  for  the  lower  Reynolds  number  runs  C  and  D,  the  region  with  the  steeper  slope 
is  between  t=2.15  and  2.32  for  Run  C3  and  t=2.02  and  2.2  for  Run  D3. 

We  emphasize  that  in  all  these  runs,  we  obtain  regions  of  different  slope  as  the  time 
evolution  of  and  6  are  examined.  Since  such  similar  trends  are  also  observed  in  Runs 
C  and  D  which  are  much  better  resolved  than  A  and  B,  we  conclude  that  the  results 
cannot  be  an  artifact  of  a  possible  resolution  problem. 

The  rate  at  which  6  decreases  is  almost  the  same  in  three  runs  (A3,  B3,  D3)  and 
slightly  slower  in  Run  C3  (Figure  25).  The  critical  times  (Tc)  is  at  a  later  time  as 
Reynolds  number  is  increased.  Similar  results  are  obtained  in  vortex  reconnection 
runs  We  make  the  following  observations:  i)  In  two  diffe;-  mt  studies  by  Brachet 

et  al.  6  evolution  with  time  seems  to  exhibit  an  exponential  decrease.  However, 

there  are  hints  that  such  a  trend  might  change  after  a  certain  time.  Brachet  et  al.  al¬ 
ways  considered  the  formation  of  a  finite  time  singularity  possible,  stating  that  ‘A  more 
singular  behaviour  due  to  subsequent  instabilities  of  the  above  structures  is  however  not 
excluded’  Such  an  interval  with  a  faster  than  exponential  rate  is  what  we  observe  in 
our  simulations,  ii)  We  find  that  S  has  a  non-zero  minimum  in  our  simulations.  This 
suggests  that  the  solution  is  uniformly  analytic.  This  is  the  case  for  the  viscous  Burgers 
equation  which  also  gives  a  6  decrease  to  a  minimum  and  a  subsequent  linear  increase 
(Figure  3  in  Reference  11). 

Next,  we  focus  on  the  spatial  location  of  the  local  maximum  vorticity.  We  find  that 
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in  all  our  runs  the  maximum  is  located  at  the  zero  planes  of  the  flattened  vortex  tubes 
and  more  specifically  in  the  “head  region”  of  the  vortices.  We  present  in  Figures  26-31 
the  magnitude  of  vorticity  contours  at  different  times  for  runs  A3  and  B2  (Recall  that  we 
have  presented  the  pictures  for  Run  A  at  earlier  times  in  Figure  13).  We  start  with  Run 
B2,  Figures  26-28.  We  observe  that  the  vortex  head  deforms  in  a  self-similar  way  until 
t= 1.937.  After  this  time,  for  example  at  t=2.0,  the  self-similarity  is  lost  in  the  sense  that 
the  vorticity  contour  pictures  can  not  be  generated  by  only  stretching  and  compressing 
the  pictiures  at  an  earlier  time.  For  Run  A3,  a  similar  picture  appears  with  a  slight 
delay  in  time.  As  seen  from  Figures  29  and  31,  vorticity  contours  seem  to  be  self-similar 
until  t=2.0.  After  this  time,  the  self-similarity  is  lost  (for  example;  the  contours  at 
t=2.062.  Figure  31).  To  quantify  this  process,  we  have  calculated  the  distance  between 
the  location  of  the  maximum  vorticity  and  the  origin  for  both  runs.  The  results  are  given 
in  Figure  32  which  displays  the  natural  logarithm  of  the  distance  as  a  function  of  time. 
We  see  in  Figure  32  that  for  both  runs,  there  is  an  interval  when  the  distance  decreases 
linearly  (in  logarithmic  scale)  which  implies  an  exponential  shrinkage  with  time.  This 
interval  is  until  t=1.939  for  Run  A3  and  t=1.9015  for  Run  B2.  But  the  curves  do  not 
continue  with  the  same  trend,  and  after  the  exponential  period  a  saturation  is  observed. 
This  is  an  interval  between  t=l. 939-2.0015  for  Run  A3  and  t=1.9015-1.939  for  Run  B2. 
To  see  the  process  clearly,  we  present  in  Figure  33  the  trajectory  connecting  the  location 
of  the  points  with  maximum  vorticity  at  different  times  (equally  spaced  in  time)  for 
both  runs.  It  is  seen  from  this  figure  that  the  approach  of  the  maxima  towards  the 
origin  halts  for  some  time  during  the  evolution.  This  smprising  “saturation”  is  similar 
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to  the  “bouncing-off”  process  of  two  vortex  rings  which  are  shot  towards  one  another. 
In  this  interaction,  there  is  a  minimum  intercentroid  distance  of  the  vortices  after  which 
the  vortices  change  direction  and  divert  their  paths.  The  saturation  process  might  also 
be  related  to  a  competition  between  the  stretching  and  dissipation  processes.  There  is 
an  interval  when  these  forces  approximately  balance  each  other  and  a  local  equilibria  is 
attained  as  in  the  case  of  viscous  Burgers  vortex.  However,  it  looks  like  the  equilibria  is 
not  stable,  the  structure  breaks  down. 

Now,  it  might  be  interesting  to  compare  the  results  of  the  time  scale  analysis  we  have 
at  different  sections  of  this  work:  For  Run  A3,  we  observe  a  faster  than  exponential 
shrinkage  (in  distance  to  origin)  for  t  >  2.0015.  The  self-similarity  is  broken  around 
this  time  (Figures  30,  31),  the  analjdicity  strip  6  has  its  fastest  decrease  rate  (between 
t=1.98-2.1),  and  a  possible  singular  time  is  predicted  to  be  at  Tc  =  2.1077.  Similarly  for 
Run  B2,  we  find  a  faster  than  exponential  shrinkage  rate  for  t  >  1.95,  the  self-similarity 
is  broken  at  this  time  (Figure  31),  £  attains  its  steepest  slope  between  t= 1.964-2.064,  and 
a  possible  singularity  time  is  predicted  to  be  at  Tc  =  2.0389.  In  spite  of  all  these  trends 
towards  a  singularity,  both  S  and  u?mai  remain  bounded  after  the  violent  deformation 
period.  We  propose  the  following  scenario:  i)  The  6  vortex  dipoles  deform  one  another, 
their  cores  get  more  and  more  flattened,  the  distances  shrink  exponentially,  ii)  There 
is  a  short  local  equilibrium  interval  when  the  viscous  forces  and  vortex  stretching  forces 
balance  each  other,  iii)  The  intermediate  state  is  likely  to  be  unstable,  the  vortices  start 
to  break  up  by  an  instability  mechanism  which  is  not  well-understood,  but  a  viscosity- 
induced  stability  is  not  discarded,  iv)  Viscous  dissipation  t2dces  over,  trends  towards 
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a  singularity  are  damped  down,  v)  In  the  runs  with  higher  Reynolds  numbers,  we  do 
not  observ  vortex  reconnection.  In  Run  C  (Re  =  500)  we  do  see  evidence  of  multiple 
reconnections,  but  further  analysis  is  necessary.  If  vortex  reconnection  does  occur,  it 
would  be  a  mechanism  against  the  formation  of  a  singularity  by  changing  the  extremely 
flattened  sheet-like  vortices  back  to  their  original  tubular  topology  and  hence  decreasing 
the  existence  of  the  high  gradient  regions. 

One  final  check  to  be  made  to  strengthen  the  above  argument  will  be  to  determine  at 
which  length  scales  the  viscous  time  scales  will  be  comparable  to  the  convective  (defor¬ 
mation)  time  scales  in  the  problem.  For  this,  we  refocus  on  Figure  32  and  particularly 
on  the  region  where  there  is  exponential  shrinkage  in  length.  For  Run  A3,  this  is  an 
interval  of  At  2;  1.939  -■  1.639  =  0.3  and  for  Run  B2,  1.9015  —  1.7015  =  0.2.  Hence, 
we  will  use  these  time  intervals  as  the  convective  time  scale  in  the  deformation  pro¬ 
cess.  The  viscous  time  scale,  on  the  other  hand,  is  given  by  l/(t/A:^).  Equating  the 
convective  time  scale  Tamv  (=0.2  for  Rim  B2  and  =0.3  for  Run  A3),  to  the  viscous  time 
scale,  we  obtain  the  wavenumber  at  which  the  dissipative  forces  have  comparable  time 
scales  as  the  convective  forces.  This  is  found  to  be  k  =  82  for  Run  A3  and  k  =  158 
for  Run  B2.  The  corresponding  length  scales  can  be  evaluated  using  L  =  2nfk  and 
one  obtains  L  =  0.0766  for  Run  A3  L  =  0.0398  Run  B2.  The  grid  size  is  given  by 
Ax  =  (7r/2)/256  =  0.0061.  Hence,  the  dissipative  time  scale  becomes  comparable  to 
convective  time  scale  at  13  gridpoints  for  Run  A3  and  7  gridpoints  for  Run  A2.  To  show 
the  resolution  clearly,  we  display  in  Figure  28  (Run  B2)  and  Figure  31  (Run  A3)  the 
vorticity  contour  sections  at  x  =  0  plane.  For  Run  B2,  we  see  that  7  gridpoint  limit 
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is  reached  somewhere  between  t  =  1.937  amd  2.0.  For  Run  A3,  the  13  point  limit  is 
reached  somewhere  between  i  =  2.0  and  2.062.  Once  these  limits  are  reached,  viscous 
dissipation  forces  will  play  a  role  against  the  singularity  formation  as  rapidly  as  the 
convective  forces  play  towards  a  singularity. 

Finally  in  this  section,  we  will  examine  the  evolution  of  enstrophy  and  skewness  in 
more  detail.  It  has  been  shown  by  Kida  and  Murakami  ^  (Figure  2)  and  confirmed  by  the 
present  work  that  the  time  at  which  the  enstrophy  attains  its  maximum  value  depends 
weakly  on  Reynolds  number.  For  higher  Reynolds  number,  the  time  of  the  enstrophy 
peak  is  later  than  that  of  the  lower  Reynolds  number,  but  the  difference  is  slight.  The 
Reynolds  number  based  on  Taylor  microscale,  R\  is  initially  674  for  Run  A  and  1685  for 
Rvm  B.  Our  results  are  in  agreement  with  the  recent  work  by  Herring  and  Kerr^  who 
observed  a  similar  trend  in  the  time  of  the  peak  enstrophy. 

The  initial  growth  of  the  enstrophy  for  all  runs  seems  to  be  at  least  exponential. 
This  can  be  seen  from  Figure  3,  17  and  19  which  approximately  give  a  linear  result 
when  the  enstrophy  evolution  is  plotted  on  a  log-lin  scale.  The  graphs  have  a  slight 
concavity  which  suggests  that  the  rate  might  be  slightly  faster  than  exponential.  Another 
interesting  feature  is  the  slight  dimple  observed  during  the  evolution.  This  is  around 
t=2.06  for  Rtm  B2  and  t=2.11  for  Rim  A3.  It  is  clearly  seen  in  Figure  19.  Thus,  we 
conclude  that  the  enstrophy  increase  rate  changes  after  the  break  up  of  the  structure. 
If  the  viscous  terms  are  ignored,  and  isotropy  is  assumed,  enstrophy  rate  can  be  related 
to  the  velocity  derivative  skevmess,  so  we  might  also  expect  a  change  in  the  skewness 
as  the  cause  of  the  enstrophy  rate  change.  This  is  indeed  the  case,  and  we  find  sharp 
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skewness  fluctuations  around  t=2.0640  for  Run  B2  and  t=2.1265  for  Run  A3.  The  time 
at  which  the  peak  is  attained  is  at  an  earlier  time  for  the  higher  Reynolds  number  run. 
Such  a  trend  is  also  observed  in  the  Taylor- Green  simulations  by  Brachet  et  al.  ®  (See 
their  Figure  10). 

It  is  a  common  practice  to  assume  isotropy  in  calculating  the  skewness,  since  the 
quantity  can  directly  be  related  to  enstrophy  for  the  inviscid  case  given  by 

fPTT  (26) 

We  have  compared  the  skewness  results  obtained  from  the  above  equation  and  Ek}uation 
9.  We  observed  major  differences  between  them  particularly  during  the  intermittent 
interval.  For  example,  for  Run  B2,  the  skewness  value  calculated  directly  is  0.52  at  t=4 
compared  to  0.00193  obtained  using  the  isotropic  relation.  Such  a  low  value  is  due  to 
the  fact  that  the  enstrophy  increase  rate  is  very  slow  around  that  time,  giving  extremely 
low  skewness  values.  Such  a  discrepancy  is  not  surprising  since  one  observes  that  close 
to  the  break  up,  there  are  two  local  regions  where  all  the  action  is  concentrated  (0 
and  x/2  pl2uies)  and  therefore  the  flow  is  strongly  anisotropic  and  the  skewness  values 
obtained  using  Equation  26  assuming  isotropy  do  not  give  realistic  values.  Similarly, 
Taylor-Green  flow  simulations  suggested  a  strong  anisotropy  before  the  break  up  (See 
Figure  8a  in  Reference  6). 

Finally,  we  make  an  observation  on  the  influence  of  the  resolution  on  the  skewness 
and  flatness.  We  have  concluded  before  that  Run  B1  is  underresolved  around  t=2.  A 
comparison  of  the  skewness  values  for  Run  B1  and  B2  shows  that  the  two  runs  approach 
asjrmptotically  to  different  values  after  the  break  up,  Run  Bl  approaches  to  0.44  and 
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Run  B2  to  0.5.  There  is  also  divergence  in  flatness  values.  Run  Bl  approaches  to  5.2 
whereas  Run  B2  tends  to  7.5  after  t=2.  Such  a  divergence  brings  questions  about  the 
calculated  turbulence  statistics  by  Kida  and  Murakami  after  the  break  down,  since  their 
maximum  resolution  is  as  much  as  our  Run  Bl  resolution.  One  other  possibility  is  that 
these  values  will  converge  to  the  same  value  at  the  later  times.  We  remind  the  reader 
that  Kida  and  Murakami’s  rims  are  longer  in  time  than  the  present  work.  To  check  this, 
we  continued  Run  B2  up  to  t=6  and  found  that  the  skewness  is  0.5295  and  the  flatness 
is  5.65  at  this  instant.  The  values  have  not  converged  yet  at  this  time. 

5  Conclusions 

The  present  work  is  an  attempt  to  show  the  similarities  and  differences  between  the 
high-symmetry  initial  condition  and  several  other  flows  in  literature  such  as  the  Taylor- 
Green  vortex  and  the  vortex  reconnection  problems.  We  And  the  following  similarities 
in  all  these  flows: 

1.  All  of  these  flows  examine  the  vortex  tubes  which  are  moving  towards  one  another. 
The  topology  is  simpler  in  the  vortex  reconnection  problems  since  there  are  only 
two  vortex  tubes  in  the  field.  In  both  Taylor-Green  and  high-symmetry  flows,  more 
than  one  dipole  exists  in  the  domain.  All  of  these  problems  result  in  the  extreme 
deformation  (flattening)  of  (at  least  one)  a  dipole,  which  is  moving  towards  the 
origin  during  the  deformation  process.  The  origin  is  a  vorticity  null  point  in  Taylor- 
Green  and  high-synunetry  flows  but  such  a  constraint  is  not  imposed  in  the  vortex 
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2.  During  the  interaction,  there  is  an  interval  during  which  local  quantities  such 
as  maximum  vorticity  and  strain  rate  and  integrated  quantities  such  as  velocity 
derivative  skewness  and  flatness  increase  sharply.  The  increase  is  so  sudden  and 
sharp  that  questions  arise  as  to  whether  some  of  these  quemtities  would  diverge  in 
finite  time  or  not.  This  interval  of  extreme  intermittency  takes  place  before  the 
peak  in  enstrophy  and  total  dissipation  is  reached. 

3.  For  viscous  simulations,  both  for  the  vortex  reconnection  problem'*’'”^,  Taylor- 
.  Green  vortex  and  the  high-symmetry  flow,  as  the  Reynolds  number  is  increased, 

the  amplitude  of  maximum  vorticity,  skewness  and  flatness  increase  and  the  peaks 
are  attained  at  an  earlier  time.  This  is  a  tendency  which  shows  hints  that  for  the 
inviscid  case,  all  these  quantities  will  blow  up. 

4.  During  this  intense  period,  all  the  flow  statistics  are  far  away  from  any  Kolmogorov- 
type  scaling  laws.  The  energy  spectra  behave  as  k~*  and  evolves  towards  k~^.  Flow 
is  extremely  anisotropic.  Skewness  and  flatness  values  are  much  higher  than  those 
obtained  in  fully  developed  turbulent  flow  simulations. 

There  are  also  differences  encountered  among  the  above-mentioned  studies.  We  summa¬ 
rize  these  as: 

1.  In  the  presence  of  viscosity,  no  singularity  formation  is  observed  in  any  of  the 
studies.  However,  both  viscous  vortex  reconnection^®’**  and  the  present  work  shows 
trends  towards  a  singularity  for  the  inviscid  case.  Such  a  trend  is  not  supported  by 
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the  Taylor- Green  simulations^*^  and  the  reasons  for  the  disagreement  are  not  clear. 
The  differences  between  the  present  work  and  the  viscous  Taylor-Green  simulations 
are;  slightly  higher  resolution  in  the  present  work  and  a  more  complicated  flow 
topology  in  the  high-symmetry  flow  compared  to  the  Taylor-Green  flow.  In  fact, 
the  violent  flow  evolution  configuration  in  high-symmetry  flow  which  is  as  a  result 
of  twelve  converging  vortices  towards  the  origin  is  not  present  in  the  Taylor-Green 
initial  conditions. 

2.  Even  though  we  observe  tendencies  towards  a  singularity,  the  presence  of  viscosity 
damps  the  process  down.  In  vortex  reconnection  problems,  vortex  line  topology 
changes  during  reconnection,  resulting  in  the  formation  of  new  tubular  vortices. 
The  flattened  remnants  are  rapidly  dissipated  by  enhanced  dissipation^^*^^.  In  the 
present  work,  the  structure  breaks  down  followed  by  a  similar  dissipation  process 
acting  against  the  singularity  formation. 

It  is  interesting  to  see  whether  the  violent  peaks  observed  in  the  simulation  of  the  high- 
symmetry  flow  can  be  avoided  by  slightly  perturbing  the  flow  field.  To  check  this,  we  have 
added  some  random  noise  (5  %)  to  the  Run  C2  at  t=2.0  and  also  at  t=2.25  (separately) 
and  continued  the  flow.  We  compared  the  local  maximum  vorticity,  skewness  and  also 
the  vorticity  isosurfaces  and  contours  (at  0  and  x/2  planes)  between  the  two  runs  with 
and  without  the  noise.  We  did  not  observe  any  significant  change  in  our  results. 

A  natural  question  is  to  discuss  what  influence  the  symmetries  have  on  our  results. 
We  tend  to  believe  that  flows  like  Taylor-Green  flow  and  high-symmetry  flow  with  all 
the  symmetry  constraints  are  unstable  and  in  the  absence  of  the  symmetries,  the  flow 
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will  escape  from  the  singularity  formation  direction.  We  remind  the  reader  that  such 
a  trend  due  to  axisymmetry  is  edso  suggested  by  Pumir  emd  Siggia*^.  It  is  not  easy  to 
check  this  since  performing  a  well-resolved  simulation  without  imposing  the  symmetries 
requires  many  more  modes  than  the  current  computers  can  provide.  But  it  has  been 
observed  that  marching  towards  or  away  from  a  singularity  strongly  depends  on  how  the 
initi2il  conditions  are  set  up^^ .  Even  the  numerical  method  chosen  might  have  influence 
on  the  results.  Finite  difference  solutions  which  do  not  make  use  of  any  symmetry  do  not 
give  any  singularity  formation^^  whereas  Fourier-Chebyshev  method  simulations  which 
impose  symmetry  constraints  intrinsically  in  the  code  obtain  a  singular  trend^'  for  the 
same  problem. 
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Table  No.I 


High-symmetry  Runs  in  Literature 


Author(s) 

Re  =  1/1/ 

Rx 

Total 

Resolution 

Independent 

Modes 

Orszag  ® 

300 

37 

32^ 

8® 

16 

Bracket  et  al.  ® 

3000,  00 

110 

256® 

42® 

84 

Bracket  ^ 

5000 

140 

864® 

144® 

288 

Kida  and  Murakami  ^ 

2000 

90 

512® 

85® 

171 

Kida  and  Murakami  ^ 

5000 

143 

512® 

85® 

171 

Present  Work 

5000 

143 

1024® 

171® 

341 
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Table  No.II 


Parameters  of  Study 


Run 

Re  =  \/v 

A# 

Resolution 

(fi?  ^/) 

A1 

2000 

0.001 

256’ 

(0,2.25) 

A2 

2000 

0.001 

512^ 

(0.785,2.25) 

A3 

2000 

0.0005 

1024^ 

(1.625,2.25) 

B1 

5000 

0.001 

512^ 

(0,6.0) 

B2 

5000 

0.0005 

10243 

(1.625,6.00) 

Cl 

500 

0.001 

2563 

(0,2.5) 

C2 

500 

0.001 

512’ 

(0.8, 2.5) 

C3 

500 

0.001 

10243 

(1.55,2.5) 

D1 

1000 

0.001 

2563 

(0,2.25) 

D2 

1000 

0.001 

5123 

(0.8,2.25) 

D3 

1000 

0.001 

10243 

(1.55,2.25) 

Table  No.III 


Memory  and  Time  Requirements 
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Resolution 


Machine 


Memory-Mwords  Cpu  secs./timestep 


256* 

Cray  2 

2.9 

2.0 

256* 

Cray  Y-MP 

2.9 

0.8 

512* 

Cray  2 

18.5 

10.0 

1024* 

Cray  2 

139.5 

80.0 

Figure  Captions 


Figure  1:  Vorticity  magnitude  isosurfaces  for  Run  B  at  t=0,  0.2  (top  row,  u  =  0.3wo), 
0.5,  0.75  (bottom  row,  Cj  —  0.49wo). 

Figure  2:  Same  as  previous  figure,  different  view. 

Figure  3;  Run  Bl,  Top:  Energy  and  enstrophy.  Bottom:  Dissipation  and  Taylor 
Reynolds  number. 

Figure  4:  4a  and  b:  Run  Bl,  energy  spectra  at  t=6.  Fig.4a:  E(k).  Fig.4b:  k^^^E{k), 
Fig.4c:  Run  B2,  k^^^E{k)  at  six  different  times  between  t=5.75  and  6.0. 

Figure  5:  Run  Bl,  Kolmogorov  wavenumber. 

Figure  6:  Symmetries  at  a  XY  plane,  z=rr/2.  Top:  u?*  and  uiy.  Bottom:  u;,  and  magni¬ 
tude  of  vorticity.  Run  A2,  t=1.25. 

Figure  7:  Vorticity  magnitude  isosurfaces  for  Run  B  at  t=1.25,  1.5,  1.75,  1.8125, 
Cf  —  1.25a)t)). 

Figure  8:  Same  as  previous  figure,  different  view. 

Figure  9:  Vorticity  magnitude  isosurfaces  for  Run  B  at  t=1.875,  1.9375  (top  row, 
u)  =  2.5u^),  2.0,  2.0625  (bottom  row,  u)  =  6.25a;b)- 
Figure  10:  Same  as  previous  figure,  different  view. 

Figure  11:  Vorticity  magnitude  isosurfaces  for  Run  B  at  t=2.125,  2.5,  3.0,  3.5  (ut  = 
6.25ci;o). 

Figure  12:  Same  as  previous  figure,  t=4.0,  6.0. 

Figure  13:  Magnitude  of  vorticity  contours  at  z=0  section  for  Run  A,  t=0.75,  1.0,  1.25, 
1.5.  A  similar  picture  exist  for  x=0  and  ysO  planes. 
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Figure  14:  Magnitude  of  vorticity  contours  at  z=itl2  section  for  Run  A,  t=1.25,  1.5, 
1.625,  ?  .75.  A  similar  picture  exist  for  x=ir/2  and  y=ir/2  planes. 

Figure  15:  Same  as  previous  figure,  at  t=1.8l25,  1.9375,  2.0,  2.0625. 

Figure  16:  Same  as  previous  figure,  at  t=  2.125,  2.5. 

Figure  17:  Evolution  of  energy  E,  enstrophy,  fi,  Kolmogorov  wavenumber  and  Taylor 
Reynolds  number  R\  for  Run  Bl  (dashed)  and  Run  B2  (solid). 

Figure  18:  Evolution  of  velocity  derivative  skewness  S,  flatness  F  and  maximum  vorticity 
^maxi  for  Run  Bl  (dashed)  and  Run  B2  (solid). 

Figure  19:  Comparison  of  two  high-resolution  runs.  Solid:  Run  B2.  Dashed:  Run  A3. 
Velocity  derivative  skewness  S,  flatness  B,  maximum  vorticity  Urnaxi  energy  E,  enstrophy 
(I  (Lin-Lin),  enstrophy  (Log-Lin) 

Figure  20:  evolution  to  determine  the  trend  towards  a  singularity,  Run  B2  (solid), 

Rtm  A3  (dashed).  Top:  Before  dissipation  takes  over.  Bottom:  After  the  peak  vorticity 
is  reached.  Fitting  for  the  top  figure  gives  Tc  =  2.039  for  Run  B2  and  Tc  =  2.108  for 
Rim  A3. 

Figure  21:  evolution  to  determine  the  trend  towards  a  singularity,  Run  D3  (solid). 

Run  C3  (dashed). 

Figure  22:  Scaling  of  Umax  and  {Tc  —  T)~^  for  Run  B2  (solid,  Tc  =  2.039)  and  for  Run 
A3  (dashed,  Tc  =  2.108). 

Figure  23:  Calculated  and  fit  energy  spectra  at  t=2.30  for  Run  C3.  Solid:  Z.5  <  k  <  340, 
Dashed:  9.5  <  ib  <  340,  Short  Dashed:  Simulations. 

Figure  24:  Energy  spectra  at  t=1.75  (pluses)  and  t=2.25  (triangles)  for  Run  B2.  Figure 
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24b:  k^E.  Figure  24c:  k^E.  The  spectra  evolves  from  k~*  to  k~^  during  the  break  up. 
Figure  25:  Evolution  of  the  width  of  the  analyticity  strip  6  for  Run  A3  (long-dashed), 
B2  (solid),  C3  (dot-dashed),  D3  (short-dashed),  (Type  2  fitting).  The  dashed  horizontal 
line  indicates  twice  the  effective  grid  size.  Log-Log  scale. 

Figure  26:  Vorticity  contours  at  z=0  (location  of  peak  vorticity)  at  t=1.75,  1.8125,  1.875, 
1.9375,  Run  B2.  The  ratio  of  the  axes  in  each  figure  is  drawn  to  scade.  But  note  that 
the  last  two  figures  have  coordinates  half  of  the  first  two  for  clarity.  The  numbers  in  the 
axes  show  the  grid  points. 

Figure  27:  Same  as  previous  figure  at  t=2.0,  2.0625. 

Figure  28:  Same  as  Figures  28  and  29  the  x  axis  stretched  to  show  the  structures  in 
more  detail. 

Figure  29:  Vorticity  contours  at  z=0  (location  of  peak  vorticity)  for  Run  A3  at  t=1.875, 
1.9375,  2.0. 

Figure  30:  Same  as  previous  figure,  at  t=2.0625,  2.125. 

Figure  31:  Same  as  Figures  31  and  32,  the  x  axis  stretched  to  show  the  structures  in 
more  detail. 

Figure  32:  Evolution  of  the  distance  from  the  origin  of  the  location  of  the  maximum 
out-of-plane  vorticity  (u;...).  Solid:  Run  B2.  Dashed:  Run  A3. 

Figure  33:  The  location  of  the  maximum  out-of-plane  vorticity  (w,)  between  t=1.7015- 
2.0640  Run  B2  (left  figure)  and  between  t=1.6390-2.2140  Run  A3  (right  figure).  Each 
marker  is  equally  spaced  in  time. 

Figure  34:  Energy  spectra  for  Run  Cl  (solid),  C2  (dot-dashed)  and  C3  (dot)  at  t=2.35. 
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Figure  35:  Evolution  of  Kolmogorov  wavenumber  fc,  for  Run  Cl  (squares),  C2  (trian¬ 
gles)  and  C3  (circles). 

Figiue  36:  Evolution  of  enstrophy  for  Run  Cl  (squares),  C2  (triangles)  and  C3  (circles). 
Figure  37:  evolution  to  determine  Tc  for  Run  Bl  (dashed)  and  Run  B2  (solid). 

Figure  38:  Evolution  of  for  Run  Cl  (squares),  C2  (triangles)  and  C3  (circles). 
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Appendix  A 


Examiiiing  a  singularity  by  computation  is  anaJogous  to  a  classical  tragic  drama;  one 
knows  from  the  stau:t  that  eventually  the  loss  of  r^lution  will  knock  the  solution  off 
the  singular  evolution  trajectory  and  doom  the  calculation. 

When  one  attempts  to  integrate  through  a  singularity,  convergence  tests  will  obvi¬ 
ously  fail,  but  how  they  fail  may  give  an  indication  of  the  type  of  singularity.  As  an 
example,  consider  the  approximation  of  a  delta  function  with  a  finite  Fourier  series.  The 
Fourier  coefficients  are  all  equal  to  a  constant  which  is  denoted  by  C.  The  delta  function 
can -be  written  as  follows: 

^(x)  =  lim  Cfj{x)  (27) 

j~*co 

where 

fjix)  =  13  exp(-*»^)-  (28) 

n=0 

and  i  =  y/^.  Since  fj{x)  is  a  geometric  series,  its  sum  can  be  found  as 

fj{x)  =  [1  -  exp(-ijJVx)]/[l  -  exp(-ix)].  (29) 

The  maximum  absolute  value  of  j  =  0, 1, ...  occurs  at  x  =  0.  In  the  limit  as  x  — ♦  0 
and  through  the  use  of  I’Hopital’s  rule,  one  finds 

ll/i(a:)||oo  =  Jim/i(a:)  =  jN.  (30) 

The  expression  above  shows  that  if  one  is  tr}dng  to  resolve  a  delta  function  with  a  finite 
Fourier  series,  each  doubling  of  the  number  of  modes  or  resolution  will  cause  doubling 
of  the  maximum  value.  Similarly,  the  maximum  value  upon  doubling  of  resolution  for 
the  representation  of  the  nth  derivative  of  a  delta  function  will  scale  as  (jW)""*"^. 
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While  the  doubling  behavior  at  the  critical  time  that  we  have  observed  in  our  simula¬ 
tions  is  a  necessary  condition  for  the  true  function  to  be  a  delta  function,  it  is  not  proved 
to  be  a  sufficient  condition.  Other  generalized  functions  may  has  a  similar  scading,  and 
a  function  that  is  continuous  but  appears  to  be  a  delta  function  in  the  resolution  range 
over  which  the  convergence  is  investigated  may  also  have  this  behavior. 
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Appendix  B:  Resolution  Checks 


In  this  section,  we  present  the  resolution  checks  performed  for  each  run  at  different 
instants.  We  systematically  checked  the  following  quantities  in  our  runs: 

1.  Spectral  convergence  for  different  resolutions.  (Energy  spectra  for  three  different 
resolutions). 

2.  High  wave-number  end  of  the  spectra,  to  see  whether  there  is  a  flip-up  in  the 
energy  in  this  region. 

3j  Comparison  of  the  maximum  resolved  wave-number  knax  with  the  Kolmogorov 
wave-number  k,. 

4.  Convergence  checks  for  total  energy,  enstrophy,  dissipation  and  the  Reynolds  Num¬ 
ber  based  on  Taylor  microscale  (for  different  resolutions). 

5.  Convergence  checks  for  local  maximum  vorticity,  skewness  and  flatness  (for  differ¬ 
ent  resolutions). 

6.  Evolution  of  6  for  different  runs. 

7.  Physical  picture  check.  Checking  the  number  of  gridpoints  across  the  vortices  in 
vorticity  contours. 

8.  Physical  picture  check.  Checking  the  vorticity  contours  and  isosurfaces  for  different 
resolutions. 
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The  Checks  1  and  2  showed  that  all  runs  are  spectrally  converged.  The  flip-up  in  the 
spectra  is  observed  to  be  the  highest  around  the  time  of  the  singularity  amd  decreases 
as  approached  to  the  time  of  peak  enstrophy.  As  an  example,  we  present  in  Figure  34 
the  energy  spectra  for  Run  C  (for  3  different  resolutions)  at  t=2.35  (close  to  the  peak 
in  maximum  vorticity). 

Check  3  showed  that  kmax  is  axound  240,  140,  89,  55  for  Rims  B,  A,  D  and  C  respec¬ 
tively  around  the  time  of  peak  local  vorticity  (Recadl  that  the  maximum  wavenumber 
for  the  highest  resolution  is  341).  These  values  increase  further  as  the  flow  approaches 
to  peak  in  enstrophy  (for  example  323  for  Run  B2).  The  check  shows  that  Run  B2  is  in 
the  limit  of  resolution  (323  vs  341)  while  Runs  A3,  C3,  D3  are  sufficienly  resolved.  We 
present  in  Figure  35  the  Kolmogorov  wave-number  for  three  different  resolutions  in  Run 
C.  We  find  that  the  three  rims  differ  by  less  than  1  %. 

Check  4  for  Run  B  was  already  mentioned  in  the  text;  we  recall  that  the  comparison 
of  the  different  resolutions  shows  that  the  energy  differs  by  less  than  1%,  enstrophy  less 
than  5%.  This  error  is  even  less  for  the  lower  Reynolds  number  runs.  The  results  for  C 
are  shown  in  Figure  36. 

Check  5  for  Runs  A  and  B  have  given  discrepancies  in  runs  of  different  resolutions 
as  mentioned  before.  We,  in  fact,  find  that  all  the  values  increase  with  the  same  rate 
up  to  a  certain  time  after  which  the  increase  further  continues  for  the  highest  resolution 
run  which  reaches  to  a  higher  peak.  While  calculating  the  Tc  values  for  these  runs,  we 
plotted  the  as  a  function  of  time  for  two  resolutions.  The  tangent  is  drawn  in  such  a 
way  that  the  line  passes  through  a  region  with  the  maximum  number  of  common  points 
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possible.  Figure  37  gives  the  plot  for  runs  B1  and  B2.  The  scaling  exponents  found  as 
such  showed  that  the  powers  for  the  two  cases  are  close  to  one  another.  For  Runs  C  and 
D,  a  better  convergence  is  attained  and  the  convergence  in  local  quantities  is  within  10 
to  20  %  respectively  (between  C2  and  C3,  D2  and  D3).  As  an  example,  we  present  in 
Figure  38  the  local  maximum  vorticity  evolution  for  Runs  Cl,  C2  and  C3. 

Check  6  showed  that  the  time  evolution  of  S  is  very  similar  in  different  runs,  the 
regions  of  different  slopes  are  observed  in  all  cases.  It  is  also  seen  that  S  values  from  run 
B3  is  very  close  to  the  double  resolution  run.  The  6  values  are  calculated  using  only  the 
1024®  resolution  runs  because  of  having  a  much  wider  exponential  high  k  region  in  these 
runs. 

Check  7  has  shown  that  the  number  of  grids  used  to  resolve  a  flattened  vortex  is 
around  7  for  Run  B2,  12  for  Run  D3  and  17  for  Run  C3. 

Check  8  showed  that  both  the  isosurfaces  amd  the  vorticity  contours  show  similar 
pictures  for  ail  runs  as  the  vortices  are  flattened  and  approach  to  0  and  7r/2  corners. 
The  amount  of  small-scale  production  however  is  much  less  in  run  C3  compared  to  the 
other  three  runs. 
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